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We study a directed coupled map lattice (CML) model in d — 2 dimensions, with two degrees of 
freedom associated with each lattice site. The two freedoms are coupled at a fraction c of lattice 
bonds acting as quenched random defects. The system is driven (by adding "energy" , say) in one of 
the degrees of freedom at the top of the lattice, and the relaxation rules depend on the local difference 
between the two variables at a lattice site. In the case of conservative dynamics, at any concentration 
of defects the system reaches a self-organized critical state with universal critical exponents close 
to the mean-field values r t = 1, r B — 2/3, and r n = 1/2, for the integrated distributions of 
avalanche durations (t), size (s), and released energy (n), respectively. The probability distributions 
follow the general scaling form P(X,L) — L~°"P{XL~ Dx ), where a « 1 is the scaling exponent 
for the distribution of avalanche lengths, X stands for t, s or n, and Dx is the (independently 
determined) fractal dimension with respect to X. The distribution of current through the system is, 
however, nonuniversal, and does not show any apparent scaling form. In the case of nonconservative 
dynamics — obtained by incomplete energy transfer at the defect bonds — the system is driven out of 
the critical state. In the scaling region close to c = the probability distributions exhibit the general 
scaling form P(X,c,L) = X~ rx V(X/£ : x(cj, XL~ Dx ), where tx = a/Dx and the corresponding 
coherence length £x(c) depends on the concentration of defect bonds c as £x(c) ~ c~ Dx . 



I. INTRODUCTION 

Critical behavior in the vicinity of a second-order phase 
transition is known to be sensitive to the presence of 
quenched random defects Quenched disorder can 

cause a change of the universality class of the critical 
behavior-in the case of weak disorder, or lead to a va- 
riety of new phenomena, as in random-field systems [Q . 
The case of strong disorder, represented for example by 
the presence of random bonds in a spin system, can lead 
to a multiplicity of metastable states and frustration, as 
in spin glasses ||. 

In contrast to conventional critical phenomena in sys- 
tems which are tuned to a critical point by varying one 
or more external parameters, much recent work H has 
examined the behavior of extended open dynamical sys- 
tems (sandpile cellular automata being a good example) 
which tend to self-organize into a metastable states with 
long-range spatial and temporal correlations. Such self- 
organized criticality (SOC) in a system which is not tuned 
to a critical point can be expected to be somewhat more 
robust and less sensitive to perturbations. 

The question of relevance of such perturbations in SOC 
is of considerable interest. A dynamic renormalization 
group study of spatially continuous models 0, which 
are believed to describe some aspects of self-organization 
found in cellular automata, shows that loss of transla- 
tional invariance due to quenched disorder is a relevant 
perturbation for SOC. Similar considerations apply to 
the role of conservation laws in the evolution rules of 
such systems. Various examples studied so far |8| reveal, 
however, that having a conservation law in the dynam- 



ics is neither sufficient nor necessary condition for the 
critical state to appear. 

We explore some aspects of these questions in the 
present work, where we introduce and study a two- degree 
of freedom (or two-color) directed coupled map lattice 
model on a 2-dimensional square lattice. In our system 
there are two variables associated with each lattice site, 
and these two freedoms are practically independent in 
the evolution, except at a random fraction c of the lat- 
tice bonds, which act as quenched random defects. This 
model (termed Model B in earlier work H) provides a 
simple example of a situation wherein quenched disor- 
der acts differently from annealed disorder |lC| . Similar 
models have also been studied in the context of signal 
transmission in a neural network p| . 

The system is driven — by adding "energy", say — to 
one of the degrees of freedom at a random site at the 
top of the lattice. Instabilities can be caused when the 
difference between the two variables at a site exceeds a 
threshold value, in which case a relaxation occurs, and 
the total energy accumulated in both states at an unsta- 
ble site is transferred to the forward neighboring sites, 
creating an avalanche. Partitioning of the energy be- 
tween the states at neighboring sites depends on the kind 
of bond (defect or normal) connecting these sites with the 
unstable neighbor (see Section II for details of the dynam- 
ics). Avalanches can be characterized through the usual 
indicators such as duration t, defined as number of steps 
that the instability progresses toward the lower bound- 
ary of the lattice, and size s — the number of lattice sites 
affected. In addition, we also monitor the total energy 
released in an avalanche (or the number of relaxations) 
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n, as well as the total number of particles leaving the 
system at the lower boundary, i.e. the outflow current 
J. 

Both conservative and nonconservative situations are 
possible by altering the energy transfer at defect bonds. 
The present numerical simulations, combined with a scal- 
ing analysis of the results help in determining the condi- 
tions under which the coupled map lattice (CML) reaches 
a SOC state, and also the universality class to which the 
critical state belongs. When the dynamics is conserva- 
tive, i.e., the total energy which is removed from one 
site appears at its neighbors, the system self-organizes 
into a critical state with universal scaling exponents. On 
the contrary, if the energy transfer is incomplete at de- 
fect bonds, our results suggest that the system is sub- 
critical, with a finite coherence length which depends on 
the defect concentration c. We determine a new scal- 
ing function for the distribution of size and duration of 
avalanches. 

The case of site defects in a critical-height sandpile 
automaton has been examined in previous work |p|,|l2|, 
where we showed that the presence of nonconserving de- 
fects can lead to a loss of SOC. Site defects (either an- 
nealed or quenched) introduce a coherence length in the 
problem, which diverges as the concentration of defects 
vanishes; the relevant exponents can be derived exactly 
|l3f , since the directed abelian sandpile with defects can 
be viewed as a random branching process. 

In Section II we introduce the model and define the 
scaling forms of various distributions. Results are given 
in Section III for conservative dynamics and in Section IV 
for the case of nonconservative dynamics, followed by a 
short summary and discussion of the results, in Section V. 



II. DYNAMICAL MODEL AND SCALING 

The coupled-map lattice studied here is patterned on 
the 2-dimensional directed abelian sandpile cellular au- 
tomaton 11 which is a simple example of an exactly 
solvable system exhibiting SOC. We associate two dy- 
namical variables (hi, h 2 ), which for convenience can be 
termed energies, with each lattice site of a 2-dimensional 
directed square lattice of linear size L. The relaxation 
process at lattice site (i,j) is determined by the actual 
values of (hi, h 2 ), as follows. If the absolute value of the 
difference between hi and h 2 exceeds a critical value d c , 
i.e., 

\h x (i,j)-h 2 (i,j)\>d c , (1) 

then the site (i,j) becomes unstable and both hi and h 2 
are reset to zero, 

hi(i,j)^0;h 2 (i,j)^0. (2) 

The lattice sites are connected by bonds which can be ei- 
ther positive or negative, and these affect the subsequent 



relaxation rules differently. (We consider the situation 
when most bonds are positive, and a fraction c of nega- 
tive bonds are distributed at random and act as defects.) 
Along positive bonds, 

hi(i + l,j±) -> hi(i + 1, j±) + (hi(i,j) + h 2 (i,]))/2, 

(3) 

and along negative bonds 
h 2 (i + l,j±) h 2 (i + l,j±) + (Xhi(i,j) + h 2 (i,j))/2. 

(4) 

j± = j ± (1 — (— 1) 4 )/2 and (i + l,j±) label the two down- 
stream neighbors of (i,j). As can be seen, for A = 1 the 
total energy hi + h 2 disappearing from site (i, j) reap- 
pears in either hi or h 2 of the forward neighbors regard- 
less of the parity of the connecting bonds — the dynamics 
is conservative. If A<1, some amount of energy is lost 
along the negative bonds and the dynamics is nonconser- 
vative. This case is considered in Section IV below. 

The system is driven by adding an energy unit at ran- 
dom sites at the input at the top (i.e., along the row 
i = 1) 

hi(l,j)^hi(l,j) + l, (5) 

and allowed to proceed following the rules embodied in 
Eqs. (|l|-^) until there are no further instabilities. This 
constitutes an avalanche. Starting from an initially ran- 
dom configuration, the system eventually reaches a SOC 
state when there are avalanches of all duration 

and size scales. These can be characterized through the 
following four quantities. 

(a) The length, I, is the total distance that an 
avalanche propagates. P(£) ~ £~( 1+Q ) is the probabil- 
ity distribution of avalanches of length I. 

(b) The size, s, is the number of sites at which re- 
laxation occurs in one avalanche. D(s) ~ s~ Ts is the 
distribution of avalanches of size s or greater. 

(c) In sandpile automata distinction can be made be- 
tween the number of particles toppled and the number of 
sites involved in an avalanche Jig ] . Similarly, here we note 
that the total amount of energy released in an avalanche, 
n, is not simply proportional to the size, s, and this in- 
troduces the distribution Q(n) ~ n~ Tn for avalanches of 
released energy > n. 

(d) The duration of an avalanche, t, is described by the 
distribution P(t) ~ t- (1+Tt \ 

From finite-size scaling arguments Jl6| these distribu- 
tions should obey the following general scaling form 

P(X,L) = L- a P(XL- Dx ), (6) 

where a is the above defined exponent of the distribution 
of lengths P(t), and X represents s, n or t. Correspond- 
ingly, Dx stands for the appropriate fractal dimension, 
which is defined via following relations: 
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where <>i denotes an average over all avalanches of se- 
lected length t (i.e., avalanches which exactly terminate 
on the £' th row). 

Note that D t is the dynamic exponent (usually de- 
noted z). In the present directed model, since avalanches 
propagate only in the forward direction, the duration is 
equivalent to the length. The dynamic exponent is thus 
exactly 1, and r t = a. This need not be the case when 
the relaxation rules are more complex |J. Furthermore, 
all the exponents are not independent, since the scaling 
relations [[| 



a = D s t s = D n t„ = n D t 



(10) 



are valid for the exponents of integrated distributions as 
defined above. 

In the present work we consider quenched disorder, i.e., 
our numerical results are averaged over several sample 
lattices, each prepared by distributing defect bonds with 
concentration c. We keep the lattice configuration fixed 
for a large number of Monte Carlo steps (which is equal to 
number of events), and consider both cases of conserva- 
tive and nonconservative transfer at defect bonds. In or- 
der to minimize effects of boundaries, we chose so-called 
free boundaries in the perpendicular direction. This is 
achieved by using a lattice of size 2L x L and initializing 
avalanches between the sites L/2 + 1 and 3L/2 at top. 
According to the above dynamic rules, only escape of en- 
ergy is possible at bottom of the lattice between sites 1 
and 2L. Results of our numerical studies are presented 
in Sections III and IV. 



III. CONSERVATIVE DYNAMICS 

For A = 1 in Eq. ^, the dynamics is conservative. 
The system reaches the SOC state, which we characterize 
through the quantities enumerated in Section II. It is suf- 
ficient to consider the disorder regime of c < 0.5 owing to 
the symmetry c — ► 1 — c, hi — > h 2 . When c = this model 
reduces to a directed abelian sandpile automaton in the 
h\ degree of freedom, and the empty state in the other 
degree of freedom, h 2 = at all sites. For nonzero c, 
the SOC state is more difficult to describe. We study the 
hystogram of total energy per site, defined as E — hi + h 2 
after a relaxation event, both for c — and few values 
of c ^ 0. In the case c = it takes nonzero values at 
the interval [0, 1], however, for c ^ much larger values 
of energies occur, although with smaller probability com- 
pared to those at E = and E = 1. Spreading of the 



distribution in the presence of defect bonds c ^ indi- 
cates that the critical state is realized via multiplicity of 
configurations of the dynamic variables hi and h 2 , which 
depends on c. 

We observe however, that the concentration of defect 
bonds c ^ does not appear to affect the critical expo- 
nents characterizing the SOC state in this model. The 
distribution of lengths of the relaxation clusters P(£), 
size D(s) and number of relaxations Q(n) for two differ- 
ent concentrations of defect bonds c = 0.1 and c = 0.5 
are shown in Fig. 1. The exponents are numerically de- 
termined to be a = 1.051±0.044, r s = 0.650±0.028, and 
t„ = 0.509 ± 0.004, independent of defect concentration 
c, suggesting that this model has universal criticality. 

In Fig. Id are shown the distributions for length and 
size of avalanches Pq(£) and Do(s) are shown for the case 
c = 0,. Slopes of these curves, corresponding to the above 
defined exponents 1 + a and 1 + t s , respectively, are es- 
timated as 1.998 ± 0.019 and 1.58 ± 0.04. For r„ we find 
the value 0.528 ± 0.022. 

In Fig. 2 we present results for the time averaged size 
<s> and relaxed energies <n> in avalanches of a fixed 
(selected) length I. According to Eqs. (0)-(||), slopes of 
these curves determine the mass-to-scale ratios (fractal 
dimensions) D s and D n , which are determined to be 
D s = 1.615 ±0.007 and D n = 1.997 ±0.008. Plotted are 
the results for c = 0.5, but it was checked that the values 
of the exponents remain concentration independent. 

In order to fully characterize the self-organized critical 
state, we have determined the various exponents and the 
corresponding scaling functions. The results of the finite- 
size scaling fit according to Eq. (||) are shown in Figs. 3 
and 4, where we used a — 1.04, D s = 1.62 and D n = 2. 

In contrast to many sandpile models, where the num- 
ber of topplings (corresponding to n in our model) is 
expected to scale with a fractal exponent jL7), we find 
that the relaxation rules (|l|)-(§) lead to rather classical 
exponent D n — 2. The distributions of the size of relax- 
ation clusters D(s,L) and the distribution of number of 
relaxations Q(n, L) satisfy the finite-size scaling form (||) 
with the above determined exponents. 

The distribution of outflow current G(J,L), which is 
the current which flows over the lower boundary of the 
system, rather than having a power-law dependence on 
J, fulfills the finite-size scaling form J16[ 



G{J,L)=L-l 3 g{JL~'t>) 



(11) 



where (3 = 2</> in the stationary state. In Fig. 5 we show 
the distribution G( J, c, L) for (a) fixed concentration of 
defect bonds c = 0.5 and various lattice sizes, and (b) for 
fixed lattice size L — 192 and various concentrations of 
defects c. The distribution G( J, c, L) for finite concentra- 
tion c of defect bonds is rather localized, as opposed to 
the case c = 0, where the distribution is broad. We find 
no apparent finite-size scaling for G{J, c, L) for nonzero 
values of c . 

We end this section with a comment on numerical val- 
ues of the critical exponents. Our results suggest that 
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the values of the exponents for avalanche duration and 
for number of relaxations are close to those in the mean- 
field SOC @Jl9|— in our notation, 1 + r t » 2, and 
l + r n w 3/2. The mean- field universality class (T^| is ob- 
tained in self-organizing sandpile models where: (i) there 
is a front-like spreading of avalanches with noninteract- 
ing sites at the front, and (ii) a global constraint exists, 
which maps the dynamic model to a critical branch- 
ing process. Although our model has somewhat more 
complex evolution rules, the exponents appear to be in 
the same universality class. (Note that in our model, 
compared to the sandpile automata, there is one more 
independent exponent, i.e., 1 + t s k, 5/3.) 

IV. THE CASE OF NONCONSERVATIVE 
DEFECTS 

When A ^ 1 in Eq. (|j) the model becomes noncon- 
servative at defects bonds. By studying the probability 
distributions of duration P(t,c,L) and size D(s,c,L) of 
avalanches for various concentrations of nonconserving 
defect bonds c, we find that SOC behavior is lost as soon 
as energy conservation is lost. 

One quantity which proves to be useful in character- 
izing the present situation is the average number of re- 
laxations <n(£)> which occur in all kinds of clusters up 
to length I (which is different from the average number 
of relaxations in clusters of selected length, which was 
considered above). In the case of conservative dynamics 
this quantity exhibits a power-law Q 

<n(£)> ~ *A*(1-t„) ) (12) 

independent of the degree of disorder. 

Numerical results are obtained by fixing A = 0.9 and 
varying c. Shown in Fig. 6 are results for the average 
number of topplings <n{£)>, for the conservative case 
(open circles) as well as for A<1 and different values of 
c, where it can be seen that the power-law is lost and 
the results depend on the concentration of defect bonds 
c (the remaining lines in Fig. 6). 

We further study the effects of concentration of de- 
fect bonds on the distributions of size and duration (or 
length) of the relaxation clusters. In Figs. 7 a and b these 
distributions are shown for L = 128 and for few values of 
c. The following general scaling form is appropriate for 
the case of lattice disorder |l2j 

P(X,c,L) = X- TX T(X/t; x (c),X/L Dx ) . (13) 

Here £x( c ) is the corresponding coherence length, which, 
as is evident in Fig. 6, varies with the concentration of 
defect bonds c. By fitting the numerical data in Fig. 
7 a and b to the scaling form ( |l3| ) we determine the 
c-dependence of £x( c )> which appears to be most sat- 
isfactorily described as £t = 1/c and £ s = c~ Ds , for the 
distributions of duration and size, respectively, r t = a, 



t s , and D s being the exponent determined in Sec. III. 
The results are shown in Fig. 7 c and d. Owing to the 
factor t/L or s/L Ds in the scaling function (|l3|), we re- 
strict the fits to rather large concentrations c (see Fig. 
7 a and b) in order to satisfy the relation ^/L < 1. It 
is interesting to note that the scaling region where the 
Eq. ([l3]) applies is restricted to a finite range of values 
of c<c* in the vicinity of the point c — 0. For instance, 
the curves corresponding to c = 0.5 in Fig. 7a, b do not 
obey the scaling form (J13|), indicating that c*<0.5. In 
the limit c = the dynamics becomes conservative and 
true SOC reappears (cf. Fig. Id.). 

Similar scaling fits were introduced earlier in Ref. 0] 
in a simple critical height model with site defects. Here 
we demonstrate that in the case of more complicated re- 
laxation rules with nonconservative bond defects, the co- 
herence lengths have the same general dependence on the 
concentration of defects, namely, £,x(c) ~ c~ Dx , where, 
in principle, the values of the exponents Dx depend on 
the dynamic model. Exact expressions for the coher- 
ence lengths in the case of the directed abelian sandpile 
model with site defects have been obtained by mapping 
the model to a random branching process [0| . 



V. SUMMARY AND DISCUSSION 

In this work we have further explored the role of de- 
fects in models of self-organized criticality. In the case 
of sandpiles with random site defects IQ , the dynamics 
is altered locally at defect sites. The situation of defects 
having a more global influence is also of interest, and we 
achieve this in the present two-color random bond model. 
Our coupled map lattice has two degrees of freedom as- 
sociated with each lattice site. Disorder is present in the 
form of quenched random bond defects, which can be 
both conserving (A = 1) or nonconserving (A ^ 1). There 
is a preferred direction of transport, which serves to sim- 
plify the dynamical evolution rules, leading to a mini- 
mal model with quenched random bonds. Preliminary 
studies [p|JTo|] have indicated that annealed and quenched 
random disorder can behave differently, and our model 
is one of the simplest examples. Similar self-organizing 
coupled map lattice models have been studied to some 
extent in the literature as models of more realistic sand- 
piles p3], abstract examples of adaptive self-organizing 
systems |^l[ , vector-state models ]22| , and as nondirected 
neural networks |Tl[| . 

With conservation, A = 1, in the limit of no disorder, 
c = 0, the model becomes an abelian CML with single 
dynamic variable (energy hi), and reaches a SOC state 
in which the distribution D{hi) is nonzero at the interval 
[0, 1], and /12 = everywhere. For nonzero concentration 
of defect bonds (c / 0), additional state hi is generated 
at fraction c of lattice sites, which affects the propagation 
of avalanches. The SOC state is again reached under the 
condition that the total energy (consisting of hi + h%) ac- 



4 



cumulated at an unstable site is completely transferred 
to its neighbors. Conservation of hi by itself appears 
to be insufficient. The critical state has more complex 
structure, evidenced by a spread in the distribution of 
the values of hi and /12, the width of the distribution 
being a strong function of c. However, the critical expo- 
nents characterizing the SOC state appear to be indepen- 
dent (within numerical uncertainty) of the concentration 
of defects c. Numerical values of these exponents suggest 
that both c = and c>0 models belong to the mean-field 
SOC universality class In particular, exponents 

for the distributions of duration (or length) of relaxation 
clusters and number of relaxations are close to the ex- 
act values 1 + a « 2 and 1 + r„ « 3/2, respectively. 
Due to the more complex dynamical rules in our CML 
model when c ^ 0, we can distinguish between r n and 
the exponent of the distribution of size of avalanches, r s , 
which are equivalent in abelian sandpile models Jl5| ; the 
present results suggest that 1 + r s « 5/3. 

By modifying the relaxation rules in a way such as 
to permit incomplete energy transfer at defect bonds, 
(i.e., /12 is not conserved while hi remains conserved), 
we study the significance of the conservation on the SOC 
state. Similar to the case of sandpile automata with site 
defects, the self-organizing CML is driven out of the crit- 
ical state (cf. Fig. 7), with the concentration of defect 
bonds c playing the role of a control parameter. The co- 
herence lengths are tuned by the concentration of defect 
bonds according to the general law £x( c ) ~ c~ Dx . In 
the case of site defects, the sandpile automaton can be 
mapped onto a random branching process, and the ex- 
act expression for the £t(c) has been derived (see Ref. 
p3[ for details). We suggest that the scaling form ( |l3| ) 
together with the expression £x(c) ~ c~ Dx for the co- 
herence lengths applies for the probability distributions 
in a wider class of subcritical self-organizing systems. 

In contrast to spatially continuous models in which dis- 
order is always a relevant perturbation [Q, our numeri- 
cal results suggest that the present coupled-map lattice 
model self-organizes into a universality class of the mean- 
field SOC which is robust to quenched random bonds, as 
long as the dynamics is conservative. However, a ques- 
tion remains as if other types of defects, such as studied 
in discrete sandpile automata in Refs. p|,p3|], combined 
with different dynamic rules may lead to continuous tun- 
ing of the universality class with a parameter in CML 
models. 
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FIG. 1. Double logarithmic plot of (a) distribution of avalanche lengths P{£) vs. £; (b) integrated distribution of avalanche 
sizes D(s) vs. s; (c) integrated distribution of number of relaxations Q(n) vs. n, for two concentration of defect bonds c = 0.1 
(dotted curves) and c = 0.5 (solid curves); (d) distributions of avalanche lengths and sizes (not integrated), Po{£) vs. I and 
Do(s) vs. s, respectively, for the case c = 0. 

FIG. 2. Double-logarithmic plot of the average values of size <s>i and number of relaxations <n> e in clusters of selected 
length £, vs. £. 

FIG. 3. Double-logarithmic plot of the integrated distribution of size of avalanches D(s,L) vs. size s for several different 
values of lattice size L, (above), and the corresponding finite-size scaling plot according to Eq. (6), (below). 

FIG. 4. Same as Fig. 3 but for the integrated distribution of number of relaxations Q(n, L) vs. n for three different values 
of lattice size L, (top), and the corresponding finite-size-scaling plot according to Eq. (6), (bottom). 

FIG. 5. Double-logarithmic plot of the distribution of outflow current G(J, c, L) vs. current J for fixed concentration of defect 
bonds c = 0.5 and various lattice sizes L, (top), and for fixed L = 192 and various concentrations c as indicated, (bottom). 

FIG. 6. Double-logarithmic plot of the average number of relaxations up to distance £ from the top of pile <n(£)> vs. I for 
conservative defects A = 1 (open circles) and for A<1 and various concentrations of nonconservative defects c as indicated. 

FIG. 7. Top: Double-logarithmic plot of the distribution of length P(£, c) vs. length £ and the integrated distribution of size 
D(s,c) vs. size s of avalanches for few concentrations of nonconservative defects (from right to left: c = 0.1, 0.2, 0.3, and 0.5) 
and fixed L — 128. Bottom: Scaling plots of the same distributions for c =0.1, 0.2, and 0.3 according to Eq. (13). 
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